function [t LTSCells IB_ADCells IB_SOCells IB_ACells Synapses] = runSim(p)
%runSim runs the core of the simulation, returning the cell arrays
%   Detailed explanation goes here

global t;
t = p.t;  %-30:.01:500;
numSFCells = p.numSFCells;
numDCells = p.numDCells;
numCells = numSFCells + numDCells;

global J_SF; % LTS input
global J_IB; % IB dendridic input
global J_a; % IB axonal input
global J_s; % IB somatic input

global LTSCells;
global IB_ADCells;

global IB_SOCells;
global IB_ACells;

global LTSIndices;
global IBIndices;

global Gates;
global Synapses;

J_SF = p.J_SF;
J_IB = p.J_IB;
J_a = p.J_a;
J_s = p.J_s;

LTSCells = initLTS(numSFCells);
[IB_ADCells IB_SOCells IB_ACells] = initIB(numDCells);

LTSIndices = 1:numSFCells;
IBIndices = numSFCells+1:numCells;

Gates = ones([numCells, numCells]);

switch lower(p.connect_mode)
    case {'detached'}
        Gates(LTSIndices, IBIndices) = 0;
        Gates(IBIndices, LTSIndices) = 0;
    case {'disconnected'}
        Gates(:,:) = 0;
    case {'ib->lts'}
        Gates(LTSIndices, IBIndices) = 0;
    case {'lts->ib'}
        Gates(IBIndices, LTSIndices) = 0;
    otherwise
end

Gates(LTSIndices, LTSIndices) = Gates(LTSIndices, LTSIndices) * 5/numSFCells;
Gates(LTSIndices, IBIndices) = Gates(LTSIndices, IBIndices) * 2/numSFCells;
Gates(IBIndices, LTSIndices) = Gates(IBIndices, LTSIndices) * .95/numDCells;
Gates(IBIndices, IBIndices) = Gates(IBIndices, IBIndices) * 0/numDCells;

Synapses = zeros([length(t), numCells, numCells]);

LTSCells = [LTSCells;zeros(length(t)-1,size(LTSCells,2),size(LTSCells,3))];
IB_ADCells = [IB_ADCells;zeros(length(t)-1,size(IB_ADCells,2),size(IB_ADCells,3))];
IB_SOCells = [IB_SOCells;zeros(length(t)-1,size(IB_SOCells,2),size(IB_SOCells,3))];
IB_ACells = [IB_ACells;zeros(length(t)-1,size(IB_ACells,2),size(IB_ACells,3))];

for i = 1:(length(t)-1)
    updateLTS(i);
    updateIB_AD(i);
    updateIB_SO(i);
    updateIB_A(i);
end
    
end